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ABSTRACT 

Aims. The present study was conducted to determine the optical extinction curve for Cerro Paranal under typical clear-sky observing 
conditions, with the purpose of providing the community with a function to be used to correct the observed spectra, with an accuracy 
of 0.01 mag airmass -1 . Additionally, this work was meant to analyze the variability of the various components, to derive the main 
atmospheric parameters, and to set a term of reference for future studies, especially in view of the construction of the Extremely Large 
Telescope on the nearby Cerro Armazones. 

Methods. The extinction curve of Paranal was obtained through low-resolution spectroscopy of 8 spectrophotometric standard stars 
observed with FORS1 mounted at the 8.2 m Very Large Telescope, covering a spectral range 3300-8000 A. A total of 600 spectra 
were collected on more than 40 nights distributed over six months, from October 2008 to March 2009. The average extinction curve 
was derived using a global fit algorithm, which allowed us to simultaneously combine all the available data. The main atmospheric 
parameters were retrieved using the LBLRTM radiative transfer code, which was also utilised to study the impact of variability of the 
main molecular bands of O2, O3, and H 2 0, and to estimate their column densities. 

Results. In general, the extinction curve of Paranal appears to conform to those derived for other astronomical sites in the Atacama 
desert, like La Silla and Cerro Tololo. However, a systematic deficit with respect to the extinction curve derived for Cerro Tololo 
before the El Chichon eruption is detected below 4000 A. We attribute this downturn to a non standard aerosol composition, probably 
revealing the presence of volcanic pollutants above the Atacama desert. An analysis of all spectroscopic extinction curves obtained 
since 1974 shows that the aerosol composition has been evolving during the last 35 years. The persistence of traces of non meteorologic 
haze suggests the effect of volcanic eruptions, like those of El Chichon and Pinatubo, lasts several decades. The usage of the standard 
CTIO and La Silla extinction curves implemented in IRAF and MIDAS produce systematic over/under-estimates of the absolute flux. 

Key words, site testing - atmospheric effects - Earth - techniques: spectroscopic 



1. Introduction 

The correction for optical atmospheric extinction is one of the 
crucial steps for achieving accurate spectrophotometry from the 
ground (see Burke et al. I2010l for a recent review). Moreover, 
the study of the extinction allows to better characterise an ob- 
serving site, enabling the detection of possible trends and ef- 
fects caused by transient events, like major volcanic eruptions. 
For this reason, many studies were conducted for a number of 
observatories around the world to derive and monitor it (Tug 
1 1 9771 Sterken & Jerzykiewics 1977; Gutierrez-Moreno, Moreno 
& Cortes [19821 [19861 King [19851 Rufener[l l86l L ockwood & 
Thompson |1986t Angione & de Vaucouleurs |1986t Krisciunas 
et al. [19871 Minni ti, Cla ria & Gomez [T9891 Krisciu nas [19901 
Pilachowski et al. [T99T1 Ster ken & Manfroid [19921 B urki et 
al. [T9951 Schuster & Parrao I200T1 Parrao & Schuster 120031 

Send offprint requests to: F. Patat 

* Based on observations made with ESO Telescopes at Paranal 
Observatory. 



Kidger et al. 120031) . In the large majority of the cases this is done 
by means of multi-colour, broad- and narrow-band photometry. 
Spectrophotometry has been used in a limited number of studies, 
mainly in connection to the calibration of spectro-photometric 
standard star s (Stone & B aldwin [19831 B aldwin & Stone [19841 
Hamuy et al. [7^921 [19941 Stritzinger et al. 120051 ). 



Although several works have been published for two astro- 
nomical sites located in the Atacama Desert, i.e. La Silla (Tug 
[19771 Ste rken & Jerzykiewics [19771 Rufen er [T9861 Sterken & 
Manfr oid [19921 Grothues & Gochermann IT9921 Burki et al. 
[19951) and Cerr o Tololo (Stone & Bald win [19831 Gutierrez et 
al. 119821 119861 Stritzinger et al. 120051) . a systematic work for 
the site of Cerro Paranal (2640 m above sea level, 24.63 de- 
grees South), hosting the ESO Very Large Telescope (VLT), was 
still missing. Broad-band UBVRI extinction coefficients are rou- 
tinely obtained at this observatory since the early phases of VLT 
operations (Patat 2003 ), in the framework of quality control and 
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Table 1. PARSEC Target Stars 
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Fig. 1. Model optical transmittance computed for Cerro Paranal 
with the LBLRTM code (see Sect.[5]l for one airmass. The main 
extinction components are identified. The scale to the right refers 
to the ozone and aerosol components only. For presentation the 
model has been convolved with a Gaussian profile (FWHM=12 
A). The curves on the top panel trace the normalised UBVRI 
Johnson-Cousins passbands. 



instrument trending] In principle, as La Silla and Cerro Tololo 
are placed at similar elevations above the sea level and latitude, 
no significant differences are expected in the extinction wave- 
length dependency, so that so far the extinction curves for these 
two observatories have been widely used to correct the optical 
spectra obtained at Paranal. However, Cerro Paranal is situated 
much closer to the Pacific Ocean (12 km), and this might turn 
into a different composition of the tropospheric aerosols, result- 
ing in a different extinction law. 

Moreover, the extinction curves currently in use for these ob- 
servatories were derived more than 15 years ago, before the two 
major eruptions of El Chichdn (1982) and Pinatubo (1991), so 
that they are most likely outdated. The available broad-band data 
for Paranal are not sufficient to address these issues. Therefore, it 
was decided to undertake the PARanal Spectral Extinction Curve 
project (hereafter PARSEC), with the manifold aim of a) ob- 
taining a spectroscopic function for the extinction correction of 
optical spectra, b) getting an estimate of the variability of the 
various extinction components in a time lapse free of major vol- 
canic eruptions affecting South America, and c) setting a term 
of reference for future studies and trend analyses at this site. 
This is particularly important in the view of the construction of 
the European Extremely Large Telescope on Cerro Armazones, 
which is located only 30 km away from Cerro Paranal. In this 
article we present the results of the PARSEC project, which was 
carried on between October 2008 and March 2009. 

Very briefly, the optical extinction can be described by three 
separate components (Hayes & Latham [19751 1: Rayleigh scatter- 
ing by air molecules, aerosol scattering, and molecular absorp- 



(J2000) (J2000) 



GD108 1 10:00:47.3 -07:33:31.2 13.56 -0.22 

GD50 2 03:48:50.1 -00:58:30.4 14.06 -0.28 

EG274 3 16:23:33.8 -39:13:47.5 11.03 -0.14 

BPM16274 4 00:50:03.2 -52:08:17.4 14.20 -0.05 

EG21 5 03:10:31.0 -68:36:02.2 11.38 +0.04 

FeigellO 6 23:19:58.4 -05:09:55.8 11.82 -0.30 

GD71 7 05:52:27.5 +15:53:16.6 13.02 -0.25 

Feige67 8 12:41:51.8 +17:31:20.5 11.81 -0.34 

(1) Oke d!990t ; (2) Hamuy et al. ( H992D ; (3) Bohlin etaT 
(4) Bohlin, Colina & Finle\ fl995l 
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1 See http : //www. eso . org/observing/dfo/quality/ 



tion (O2, O3 and H2O). The contribution of the various com- 
ponents is illustrated in Fig. [TJ where we present the transmit- 
tance computed for Paranal using the LBLRTM code (Clough et 
al. [2005). While the Rayleigh scattering and aerosols act at all 
wavelengths, the effect of ozone in the optical is limited to the 
so-called Chappuis bands (Chappuis 1880) between 5000 and 
7000 A, and the Huggins bands (Huggins, [18901 below 3400 
A. The molecular bands of O2 and H2O are relevant only above 
6500 A. 



2. Observations 

2. 1 . Instrumental setup 

All observations were carried out using the FOcal Reducer/low- 
dispersion Spectrograph (hereafter FORS1), mounted at the 
Cassegrain focus of the ESO-Kueyen 8.2 m telescope 
(Appenzeller et al. |1998l Szeifert et al. l2007l) . At the time of the 
observations discussed in this paper, FORS 1 was equipped with 
a mosaic of two blue-optimised, 2048x4096 15/mi pixel (px) 
E2V CCDs. The spectra were obtained with the low-resolution 
G300V grism coupled to a long, 5.0 arcsec wide slit, giving a 
useful spectral range 3100-9330 A, and a dispersion of ~3.3 A 
px~'. The spatial scale along the slit is 0.252 arcsec px -1 . The 
G300V grism shows a rather pronounced second order contam- 
ination, starting at about 6100 A (Szeifert et al. 2007). To cover 
the widest possible wavelength range we run the observations 
using two setups, one with no filter (blue setting) and one with 
the order sorting filter GG435 (red setting). The latter provides 
a wavelength range free of second order contamination 4400- 
8100 A. 

The E2V detectors suffer from a marked fringing above 6500 
A (Szeifert et al. [2007). This becomes extremely strong for 
A >8000 A, reaching a peak-to-peak amplitude of ~20%, and 
it is not possible to remove it from the data. Therefore, the spec- 
tral region above this wavelength is practically unusable for our 
purposes. 



2.2. Target selection and observations 

In principle, there is no need to observe spectrophotometric stan- 
dards for deriving the extinction curve, as any relatively bright 
star could be used for this purpose. However, observing stan- 
dard stars has the advantage that the same data can be used also 
for the calibration plan purposes, hence mitigating the impact 
on normal science operations. For this reason, the programme 
stars were chosen among those included in the FORS 1 calibra- 
tion plan. To limit the strength of the Balmer photospheric ab- 
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Table 2. Log of the observations 

GRISM 300V (blue setting) 

Star Number Time Range Nights Airmass distribution 
Id. of spectra (MJD-54000) 1.0-1.2 1.2-1.5 1.5-2.0 >Z0 
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sorption lines (which tend to hamper the derivation of the ex- 
tinction curve, especially in the blue domain), preference was 
given to hot, blue objects. The selected target stars are listed 
in Table [T] which also presents their main properties. Exposure 
times ranged from five seconds to a minute for the faintest tar- 
gets (i.e. GD50, BPM16274). Since the maximum shutter timing 
error across the whole field of view of FORS 1 is ~5 ms (Patat & 
Romaniello 2005 ), the photometric error associated to the expo- 
sure time uncertainty is 0.1% or better. 

The spectroscopic data were collected on a time range span- 
ning about 6 months, between October 4, 2008 and March 31, 
2009. The target stars were observed randomly, mainly during 
morning twilight. On a few occasions the same star was observed 
several times during the same night, covering a wide range in air- 
mass (hereafter indicated as X). In total, more than 300 spectra 
were obtained for each of the two setups, with X ranging from 
1.0 to 2.6. The exact number of data points, time and airmass 
ranges are shown in Table |2j together with the airmass distri- 
bution for each of the eight programme stars. The observations 
were carried out under photometric or clear conditions, which 
were judged at the telescope based on the zero-points delivered 
by the available imaging instruments (transparency variations 
<10% across the whole night in the V passband). The quality of 
the nights was also checked a posteriori, using the data provided 
by the Line of Sight Sky Absorption Monitor (LOS SAM), which 
is part of the Differential Image Motion Monitor (DIMM) in- 
stalled on Cerro Paranal (Sandrock et al. 2000). For this purpose, 
we have retrieved the DIMM archival data for all relevant nights, 
and examined the RMS fluctuation of the flux of the star used by 
the instrument to derive the prevailing seeing When this fluc- 
tuation exceeded 2% (or no LOSSAM data were available), the 
night was classified as non suitable, and the corresponding data 



2 The fluctuations, measured at 5000 A, are evaluated over one 
minute time intervals. 



Table 3. FWHM image quality distribution in PARSEC spectra. 
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rejectee^] Table [2] lists only data that passed this selection (600 
out of 672 initial, non-saturated spectra). 

All spectra were obtained with the slit oriented along the 
N-S direction, which is strictly optimal only for observations 
close to the meridian. The misalignment between the slit and 
the parallactic angle at large hour angles (i.e. for targets at high 
airmass), can potentially lead to slit losses due to the differen- 
tial atmospheric refraction (Filippenko 11982b . and the seeing 
increase at larger zenith distances (Roddier ll981| ). As a conse- 
quence, the estimated extinction coefficient might be systemati- 
cally overestimated. To minimise this effect we used a 5 arcsec 
wide slit for all observations; this, coupled to the typical seeing 
attained during our observations (see Table [3), ensures that dif- 
ferential light losses are negligible for all data taken at X <2. 
Additionally, FORS1 is equipped with a Linear Atmospheric 
Dispersion Corrector (LADC), which is capable of maintain- 
ing the intrinsic image quality down to airmass ~1.5 (Avila, 
Rupprecht & Becker [T997I ). At larger airmasses the LADC only 
partially reduces the effect of the atmosphere. 

The image quality values (FWHM) reported in Table [3] were 
deduced directly from the data, analysing the profiles perpen- 
dicular to the dispersion direction at different wavelengths. The 
Table presents minimum (s mi „), maximum (s, mu -), median (s me d) 
and 95-th percentile (S95) of the observed image quality distri- 
bution. The image quality turns out to be better than ~1.5 arcsec 
in 50% of the cases. While we did not attempt to account for the 
effects of atmospheric refraction, we have applied a correction 
for the slit losses caused by seeing. The method is described in 
Appendix [A] 

3. Data reduction 

The data were processed using the FORS Pipeline (Izzo, de 
Bilbao & Larsen 2009). The reduction steps include de-bias, 
flat-field correction, and 2D wavelength calibration. The spec- 
tra were then extracted non-interactively using the apall task 
in IRAfj^jafter optimizing the extraction parameters. Because of 
the wide (5 arcsec) slit used to minimise flux losses, the uncer- 
tainty in the target positioning within the spectrograph entrance 
window is expected to produce significant shifts in the wave- 
length solution. For this reason, for each standard star we se- 
lected a low-airmass template spectrum (drawn from the data 
sample), to which we applied a rigid shift in order to match O2 
and H2O atmospheric absorption bands^] Then, the shift to be 
applied to each input spectrum was computed via cross corre- 
lation to the corresponding template spectrum. For doing this 
we first subtracted the smooth stellar continuum estimated by 

3 The experience accumulated in Paranal shows that RMS fluctua- 
tions larger than 2% on the timescale of minutes are most likely associ- 
ated to the presence of clouds. 

4 IRAF is distributed by the National Optical Astronomy 
Observatories, which are operated by the Association of Universities 
for Research in Astronomy, under contract with the National Science 
Foundation. 

5 The exposures were too shallow to show night sky emission lines. 
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the continuum task in IRAF. Wavelength shifts exceeding 10 A 
(~3 px) were recorded during this process. 

The RMS error of the wavelength solution is of the order of 
1A. Since the bluest line used by the FORS1 pipeline is Hei 
3888A, systematic deviations at shorter wavelengths could in 
principle occur. To check this possibility, we inspected the cali- 
bration frames after applying to them the same wavelength solu- 
tion used for the programme stars. Thanks to the enhanced blue 
sensitivity of the EEV detector, several lines (identified as Hg i 
3650 A, Cdi 3610 A and Cdi 3466, 3468A), are clearly de- 
tectable. The measured wavelengths of these features turn out 
to be within ~3A from their laboratory values, corresponding 
to ~l/4 of the typical FWHM resolution. This indicates that the 
wavelength solution is accurate to within a few A down to ~3400 
A. To quantify the effect of wavelength calibration inaccuracies 
of this order, we derived the extinction curve applying a wave- 
length offset to the data. For an offset of +3.3 A (one pixel), 
the variation is below +0.001 mag airmass -1 in the red, while it 
reaches +0.003 mag airmass -1 at the blue edge (~3400 A). 

Since the adopted slit width is much wider than the typi- 
cal seeing characterizing our data set, the actual resolution in 
the spectra depends on the prevailing seeing conditions. While 
this has no appreciable influence on the smooth stellar contin- 
uum, Balmer lines are expected to be affected, especially in their 
cores. In principle, this can be mitigated by matching the reso- 
lutions of the input spectra by a Gaussian profile convolution. 
However, in sight of the targeted resolution of the extinction 
curve, which is several tens of A (see Sect. and the fact that 
different stars have distinct line profiles and depths anyway, we 
have not attempted to correct for this effect. As a consequence, 
weak spurious features are expected at the wavelengths corre- 
sponding to the strongest hydrogen lines, such as Ha and H/3. 

During the PARSEC campaign, no instrument intervention 
or mirror re-aluminization has occurred, guaranteeing the homo- 
geneity of the data. However, it is known that the main mirrors 
of the VLT are subject to a reflectivity loss, mainly due to alu- 
minum oxidation. The losses are as large as ~13% per year in the 
U band, they decrease at redder wavelengths, and reach ~5% per 
year in the / band (Patat 2003). Since the data presented in this 
paper span over six months, the implied efficiency variation is 
significantly higher than the measurement errors, and this is ex- 
pected to increase the noise. To account for this effect, at least 
to a first order, the measured instrumental magnitudes were cor- 
rected using the time elapsed from the previous aluminization 
(May 17, 2008; MJD=54603.5), and a linear interp olation of the 
values deduced from broad band photometry (Patat 2003 ) to the 
required wavelengths. The application of this correction reduces 
the RMS scatter from the best fit solution by ~25%. 

Finally, airmasses were derived from the unrefracted zenith 
distances by means of the classical formula proposed by Hardie 
( 1 1962b . We note that since ~83% of the spectra were obtained at 
X <2, airmasses computed using other formulations discussed 
in the literature are always to within ~0.001. 

4. Derivation of the extinction curve 

The first step in the derivation of the extinction curve is the cal- 
culation of instrumental magnitude within predefined spectral 
bins. This is defined as follows: 

pA+AA/2 

m(A) = -2.5 log s(A) dA 

JA-AA/2 




Q U I I ill I I L_l I I I I I I I J I I I I I I L 



4000 5000 6000 7000 8000 
Wavelength [A] 

Fig. 2. Best fit extinction coefficient for the blue (circles) and 
red (triangles) settings. The errorbars indicate the 3-sigma un- 
certainty of the best fit solution. The upper and lower grey curves 
trace the 95% variation level of single measurements. The solid 
curve is an LBLRTM simulation for Paranal (see Sect. BJ, while 
the vertical segments at the bottom mark the positions of the first 
7 Balmer lines. 

where s(A) is the observed spectrum (in electrons s~') and AA is 
the spectral bin size. The fluxes within the bins were computed 
by numerical integration of the extracted, wavelength-calibrated 
spectra. In the following we will always consider bins of A/l=50 
A (15 px). This choice is motivated by three reasons: a) the re- 
sulting resolution is sufficient to follow in great detail the be- 
havior of the continuum extinction; b) since the bin width is 
significantly larger than the typical spectral resolution (~15 A 
FWHM), the values measured within adjacent bins should be 
completely uncorrelated, and c) the resulting nominal signal- 
to-noise ratio per resolution element ranges from 100 to 1200 
across the whole wavelength range covered by the observations 
(95% of the bins have a signal-to-noise ratio larger than 500). 
Since below ~7000 A all spectra are dominated by the object's 
shot-noise, the precision of the instrumental magnitudes is of 
the order of 1 % or better. Above this wavelength the dominant 
source of uncertainty is fringing, which reaches peak-to-peak 
amplitudes of ~20%. This limits the usable range to wavelengths 
shorter than 8100 A. The instrumental magnitudes were com- 
puted within adjacent bins up to 6800 A. To mitigate the effect 
of fringing, at larger wavelengths we have selected three bins 
(having widths of 160, 200, and 200 A), centred at 7060, 7450, 
and 7940 A respectively, in order to avoid the O? and H2O bands 
(Fig.Q. 

Once the instrumental magnitudes are corrected for the effi- 
ciency degradation (see previous section) and slit losses due to 
bad seeing (see Appendix [A}, they can be used to finally derive 
the extinction curve. 

Among all possible algorithms for combining data obtained 
for different programme stars, we opted for the global solution 
proposed by Hayes & Schmidtke (1987). This method is appli- 
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cable to spectro-photometry and narrow-band photometry, i.e. 
when the passbands are nearly monochromatic, as is our case. 
Under these circumstances there are no colour-terms to be taken 
into account, and each wavelength bin can be treated indepen- 
dently from the others. One issue with this method is that, by 
construction, one cannot also solve for the zero-point of the lin- 
ear relation at the same time that the slope (i.e. the extinction 
coefficient) is being determined. However, for our purposes this 
is not a problem, since we are only interested in the extinction 
term. 

The result is presented in Fig. [2] As the values derived from 
the two settings in the intersection region (4400-6075 A) are 
fully consistent within the estimated uncertainties, we averaged 
them. Also, we replaced the values corresponding to the strong 
Balmer lines with a linear interpolation between the adjacent 
bins. In Fig. [2] we also included the 95% confidence level de- 
duced from the distribution of k(A) derived from each single ob- 
servation (the extinction variations are discussed in more detail 
in Sect. [6]). 

5. Comparison with atmospheric model 

To derive the basic physical parameters related to the extinction 
curve of Paranal we used the Line By Line Radiation Transfer 
Model (LBLRTM; Clough et al. |20U3. This codffi based on 
the HITRAN database (Rothman et al. 120091 ). has been vali- 
dated against real spectra from the UV to the sub-millimeter, 
and is widely used for the retrieval of atmospheric constituents. 
LBLRTM solves the radiative transfer using an input atmo- 
spheric profile, which contains the height profiles for tempera- 
ture, pressure, and chemical composition. The code also includes 
the treatment of continuum scattering, and has an internal model 
for tropospheric aerosols (based on LOWTRAN). 

For the atmospheric profiles we adopted a standard equato- 
rial profile as a basis for all calculations. To make the simula- 
tions more realistic we then replaced the standard profiles of 
pressure, temperature, and water content with an average ver- 
tical profile derived from the Global Data Assimilation System^] 
(GDAS). The average profile, provided to us by W. Kausch and 
M. Barden (Kausch & Barden 2010, private communication), 
was obtained averaging GDAS data over the last 6 years for the 
Paranal site. The corresponding amount of precipitable water va- 
por (PWV) is 2.0 mm, while the ozone column density is 238.8 
Dobson Units (DU^j Vacuum wavelengths were converted into 
air wavelengths using the relation derived by Morton ( 1991 ). 

Following Hayes & Latham (119751 1. rather than assuming 
an aerosol model, we derived it from the data. For doing this 
we disabled the aerosol calculation in LBLRTM, and we com- 
puted the aerosol term k aa (A) as the difference between the data 
and the resulting model. The outcome is displayed in Fig. [3] As 
first suggested by Angstrom (119291 [19641, the aerosol extinc- 
tion can be described as k ae! (A) - k A~ a , where k(, and a de- 
pend on the column density and size distribution of the aerosol 
mixture. This analytical formulation has been used in a num- 
ber of extinction studies (Hayes & Latham 1975; Angione & de 
Vauco ureurs |19861 Gutier rez-Moreno et al.|19821U986i R ufener 
[19861 Minniti et al. 1989; Burki et al. [T995l Mohan et al. [T999l 
Schuster & Parrao 2001 ), and the typical value for a is -1.4±0.2 

6 See http : //r tweb . aer . com/lblrtm . html 

7 GDAS is maintained by the Air Resources Laboratory 
(http : //www . arl . noaa . gov/ 1. 

5 One DU corresponds to an ozone column density of 2.69xl0 16 
molecules cirT 2 . 
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Fig. 3. Aerosol extinction derived from the observed data and 
the LBLRTM model. The long-dashed line is a best fit model for 
^aer(^). with £o=0.014 and a--\A2 (see text). Overplotted are 
also the Rayleigh (short-dashed), ozone (solid thin curve), and 
clean tropospheric aerosol (thick solid) components computed 
by LBLRTM (see text). 
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Fig. 4. Deviations of the derived extinction curve from the 
LBLRTM simulation for Cerro Paranal (observed minus model). 
Errorbars are at the 1-sigma level. The deviation seen at about 
6300 A corresponds to the weak O2 absorption peaking at ~6280 
A. The dotted-dashed lines indicate the 3-sigma deviation from 
the best fit solution in the range 4000-6800 A (-0.007 mag 
airmass -1 ). The vertical segments at the bottom mark the po- 
sitions of the first 7 Balmer lines. 
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(with A expressed in /urn; see for instance Burki et al. 1995 ). This 
value is common for the so-called meteorological haze (i.e. free 
of volcanic pollutants) which characterises photometrically sta- 
ble nights. 

A reasonable fit to our data above 4000 A is achieved for 
k =0.01 3+0.002 mag airmass -1 and a = -1.38+0.06. This re- 
sult is fully in line with those found for other astronomical sites 
in the Atacama desert, like Cerro Tololo (Gutierrez et al. |1982| l, 
and La Silla (Rufener [19861 Burki et al. [19951 However, the 
down-turn visible below 4000 A is probably indicating a de- 
parture from a simple power-law model. We will get back to 
this issue while discussing the comparison with other observa- 
tories in the Atacama desert (Sect. TJ). Although the properties 
of aerosols and the local conditions are difficult to model (Hess, 
Koepke & Schult [l998l . for the sake of completeness we plotted 
the LBLRTM clean tropospheric aerosol model in Fig. [3] (solid 
thick curve). In order to reach a reasonable fit we had to reduce 
by 25% the LBLRTM default aerosol column density in the first 
10 km. Even larger reductions are necessary if one uses the ma- 
rine or desert aerosol types implemented in the code. 

To obtain information on the aerosol characteristics, we have 
used the Optical Properties of Aerosol and Clouds (OPAC) soft- 
ware package (Hess et al. 1998). A best fit to the data above 
4000A indicates that the mixture is mostly composed by water 
soluble particles, with very low contamination by insoluble ma- 
terial, minerals and sea salt. 

The final model for Paranal, obtained adding the best fit 
aerosol term to the LBLRTM results for the Rayleigh and ozone 
contributions, is plotted in Fig. [2] where it is compared to the 
data. The model reproduces very well the overall shape, and it 
also accurately follows the bumps in the 5800-6000 A region re- 
lated to the ozone component and O2 (see Fig. [TJ. The quality of 
the match is better seen in Fig. [4] where we plot the residuals: the 
relative deviations are less than 0.01 mag airmass -1 for the vast 
majority of the data. The largest discrepancies are seen redwards 
of 7000 A, and bluewards of -4000 A. While those in the red re- 
gion are most likely explained by fringing, the deviations seen 
in the blue are more difficult to explain. Although most of the 
data points are formally consistent with a null deviation within 
a 3-sigma level, there is a clear trend which is suggestive of a 
systematic effect (see below, and the discussion in Sect. |7]i. 

As far as the ozone content is concerned, we notice that the 
good fit achieved in the region of interest (5000 to 7000 A) in- 
dicates this is, on average, as predicted by LBLRTM for the site 
location, i.e. -240 DU. 



6. Extinction variability 

An example of the extinction time evolution is presented in the 
upper panel of Fig. ||]for /i=6000 A. Although there might be an 
hint of a mild modulation in the observed values (see the large 
dots in Fig. [5] and the discussion in Sect. IT), the limited time 
coverage of our sample (~6 months) does not allow us to firmly 
derive possible seasonal trends. Nonetheless, the available data 
are sufficient to study the extinction fluctuations under what can 
be considered as typical clear-sky conditions on Paranal. 



6. 1 . Variability of continuum extinction 

In general, the distribution of the extinction coefficient appears 
to be slightly skewed towards larger values, as found by several 
authors (see for instance Reimann, Ossenkopf & Beyersdorfer 
[19921 Kidger et al. 2003 ; Parrao & Schuster l2003l . This is illus- 
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Fig. 5. Upper panel: time evolution of the extinction at 5500 A. 
The horizontal dashed line marks the best fit value. The large 
dots are the average values computed within the time intervals 
indicated by the horizontal bars. The vertical bars are the as- 
sociated RMS deviations The time intervals include the same 
number of measurements (61), and the points are placed at the 
average time within each bin. Lower panel: distribution of the 
extinction coefficient at 4000 A (left) and 6800 A (right). 



trated in the lower panel of Fig. [5] where we present the distri- 
butions for two different wavelengths (4000 and 6800 A). The 
distribution appears to be broader in the blue than in the red. 
The peak-to-peak variations at the 95% level show a minimum 
value of 0.04 mag airmass -1 at about 6000 A, and grow to ~0.07 
mag airmass -1 at 4000 A. The inter-quartile range is ~0.01 mag 
airmass -1 for A >4000 A, while it reaches 0.04 mag airmass -1 
at the blue edge of the spectral range. 

An important aspect is that the errorbars at various wave- 
lengths are much larger than the scatter shown by the data points 
along the extinction curve (see Fig. [3J. The explanation is that 
the estimated errorbars are dominated by overall variations of 
the extinction rather than by random errors. In order to show this 
behavior in a more quantitative way, we studied the correlation 
between the extinction measured within the different wavelength 
bins. Since, for any given spectrum, k(A) is measured exactly un- 
der the same conditions, this allows us to characterise the rela- 
tion between the extinction variations across the whole spectral 
range covered by each of our setups separately. The result of this 
operation for the blue setting is presented in Fig. |6| which has 
been calculated using the Pearson coefficient r„ as a figure of 
merit for the correlation. 

Two important facts emerge from this analysis. The first is 
that the correlation is strong (r vv >0.8) above 4000 A, implying 
that the various portions of the extinction curve tend to vary in 
unison, as expected in the case of changes in the aerosol con- 
tent/composition. The second is that the extinction at A <4000 A 
appears to be less tightly correlated (r xy <0.7) to that measured 
at larger wavelengths. We believe this is only partially due to 
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Fig. 6. Extinction correlation map for the blue setting. The con- 
tours trace the iso-correlation levels. For presentation the map 
has been smoothed with a Gaussian filter (cr=75A in both di- 
rections). The vertical and horizontal strips correspond to the 
Balmer lines. The decrease of correlation seen along the diago- 
nal between 3700 and 4000 A is due to the smoothing. 

the larger random errors, and it is caused by an independent be- 
havior of the extinction curve below 4000 A (see the discussion 
in Sect. [7J. Although the correlation is always very significant 
above 4000 A, it becomes weaker for A >5700 A, in the sense 
that the extinction in the red is less correlated with that measured 
below 5500 A. One possible explanation for this is a combined 
effect of the independent variations in the aerosol (which affects 
all wavelengths), and in the ozone content (which mostly im- 
pacts the region between 5000 and 7000 A). 

6.1.1. Rayleigh scattering 

During the nights used for the PARSEC project, the atmospheric 
pressure (retrieved from the Ambient Monitor, Sandrock et al. 
2000) showed a peak-to-peak excursion of 6.5 mb around the 
average value of 743.6 mb (RMS variation is 1.1 mb). Since 
the relative variation of the Rayleigh term is proportional to 
the relative variation of atmospheric pressure (see for instance 
Hansen & Travis 119741) . this is expected to be constant to 
within 0.15% (RMS), with peak-to-peak fluctuations of less than 
1%. Therefore, larger extinction variations can be attributed to 
changes in the aerosol content and, to a smaller extent and only 
within the wavelength range 5000-7000 A, to variations in the 
stratospheric ozone column density (see Sect. 6.2.1 1. 



6.1.2. Aerosol 

To estimate the variability of the aerosol component, we have 
computed the residual extinction at 4000 A (where the ozone 
contribution can be neglected), after subtracting the LBLRTM 
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Fig. 7. Upper panel: distribution of aerosol extinction at 4000 
A. Lower panel: distribution of ozone extinction at 6000 A. The 
upper scale indicates the corresponding ozone column density in 
DU (see text). The vertical dashed lines indicate the median, first 
and third quartile of the distributions. 



Rayleigh term to the blue setting data. The resulting distribu- 
tion of the aerosol extinction, which we indicate as £ aer (4000), 
is presented in Fig. 7] (upper panel). The median value is 0.045 
mag airmass -1 , and the semi-interquartile range is 0.009 mag 
airmass -1 . Only in ~1% of the cases £ aer is smaller than 0.01 
mag airmass -1 , while it can reach values larger than 0.10 mag 



airmass (~2.5%). 



6.2. Variability of molecular bands 

In the wavelength range covered by our observations, the most 
relevant absorption bands are from molecular oxygen (O2 and 
O3) and water. While the Chappuis bands of O3 are important be- 
tween 5000 and 7000 A, O2 and water mostly affect wavelengths 
larger than 6500 A (see Fig. [I}. In the following we briefly dis- 
cuss the effect of their variability. 



6.2.1. Ozone bands 

Ozone is the main responsible for the cutoff in the atmospheric 
transmittance below 3400 A, where the O3 component quickly 
surpasses the Rayleigh scattering (the ozone term reaches ~2.5 
mag airmass -1 at 3000 A). Since our data do no cover this wave- 
length region, we estimated the ozone column density measur- 
ing the extinction ko, corresponding to the peak of the Chappuis 
bands at ~6000 A (see Fig. [3}. This was achieved removing 
the Rayleigh contribution (which was assumed to be constant 
in time), and an estimate of £ aer at 6000 A. Since the aerosol 
is strongly variable, for each data set we h ave extrapolated the 
value measured at 4000 A (cf. Sect. 6. 1 .2 1, conservatively as- 
suming that k. ier follows the usual power law with a=-1.38, 
so that £ aer (6000) 



0.6 fe aer (4000). The distribution of k 0l is 
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Fig. 8. Ozone column density for the Paranal site measured by 
the Ozone Monitor Instrument on board of AURA. The two ver- 
tical dashed lines indicate the time range covered by PARSEC 
data. The circles mark the dates when PARSEC data were ob- 
tained. 

shown in Fig. [7] (lower panel). The average value is 0.039 mag 
airmass -1 (0.002 mag airmass -1 RMS), which is slightly larger 
than the LBLRTM prediction for Paranal (0.036 mag airmass -1 ). 
The peak-to-peak variation is about 0.02 mag airmass -1 , and this 
would be sufficient to produce an enhanced extinction variabil- 
ity, at least at the Chappuis bands peak (~6000 A). However, a 
significant fraction of the variance seen in the ozone extinction is 
most likely due to measurement errors, and uncertainties in the 
estimate of the underlying aerosol contribution. 

To retrieve the implied ozone column density N(Os), we 
have run a series of LBLRTM simulations scaling the O3 profile 
to obtain a total N(03) between 60 and 480 DU. In this range 
&o 3 (6000) scales linearly with the column density, and a best fit 
to the model data gives the following relation: 

/to 3 (6000)^ 1.51 x 10 -4 N(0 3 ) 

where fco,(6000) is expressed in mag airmass" 1 , and N(03) is 
expressed in DU. From this we derive an average O3 column 
density of 258 DU (the RMS deviation is 14 DU). 

To cross-check this result, we run an independent analysis 
of the ozone variability using the satellite data provided by the 
Ozone Monitor Instrument (OMI) on board of NASA AURA0 
OMI derives the daytime ozone column density by comparing 
the amount of back-scattered solar radiation in the UV and in the 
optical. The OMI O3 column density over Paranal collected be- 
tween 2007 and 2009 clearly shows a seasonal trend (see Fig. 8j, 
with maxima attained around August, September and October 
(~300 DU), and minima reached in February, March and April 
(~240 DU). The average value during the PARSEC campaign 
was about 260 DU, which is slightly larger than the value given 
by the LBLRTM model (~240 DU; see previous section). A 



Fig. 9. Example spectrum of Feige 110 obtained with FORS 1 on 
March 13, 2006 (grism G300I+OG590; airmass X=1.03). The 
signal-to-noise ratio per pixel on the continuum is ~200. The 
main molecular absorption bands are marked. Superimposed are 
two LBLRTM simulations with PWV=2mm (dotted curve), and 
PWV=6mm (solid curve). In both cases the O2 column density 
is 3.3xl0 24 cm~ 2 . Note that at the time this spectrum was ob- 
tained, FORS1 was equipped with a Tektronik detector much 
less affected by fringing. 



closer inspection to the OMI data shows that the ozone content 
steadily decreased during the time covered by our observations. 
The peak-to-peak variation is about 25%. This turns into a max- 
imum variation of ~0.01 mag airmass -1 at ~6000 A during the 
time covered by our observations. An inspection of the time evo- 
lution of feo,(6000) shows no traces of such a steady decrease, 
and the fluctuations appear to be dominated by short timescale 
variations, partially attributable to random errors. 

6.2.2. Oxygen bands 

Molecular oxygen shows two main O2 vibrational absorption 
bands centred at 6870 A and 7605 A, usually indicated as B and 
A bands, respectively (Fig. [9j. Their typical equivalent widths 
(EW) are ~6A and ~28A, which make them easily detectable 
in low resolution spectra. To quantify the variability of O2 col- 
umn density, and its effect on the extinction, we have measured 
the EW of the B band in the red setting spectra. The A band 
is severely affected by fringing, and so no very accurate mea- 
surements were possible. However, the integrated strengths of 
the two bands are well correlated, as demonstrated by a series of 
LBLRTM simulations run for different values of the column den- 
sity N(02). Additionally, they both follow a linear dependency 
on airmass, down to X-2.5. 

The measurements clearly show the EW airmass depen- 
dency, which is well reproduced by the following best fit rela- 
tion: 



9 Data can be obtained from toms . gsfc . nasa . gov . 



EW 6m = 5.47 ± 0.04 + (2.56 ± 0.06) (X - 1) 



F. Patat et al.: Optical atmospheric extinction over Cerro Paranal 



9 



2007 



2008 



2009 



2010 



15 - 



I I I I I I I I I I I I I I I I I I I I I I I I I I I II I I I I I I I I I I I 




•a 



JiSL 

i i i_ 



500 1000 
MJD-54000 



Fig. 10. Precipitable water vapor measured in Paranal between 
2007 and 2010 (Smette et al. |2007T >. The large dots and error 
bars indicate the monthly averages and the associated RMS de- 
viations, respectively. 



where EWeno is expressed in A. With the aid of this relationship 
one can correct the observed values to zenith, and derive the col- 
umn density using a standard curve of growth procedure. For this 
purpose we computed a number of LBLRTM models varying 
N(0 2 ) between 8.4xl0 23 cirT 2 and 6.7xl0 24 cnr 2 . Subsequently 
we measured the EW of the B band on the output spectra, after 
convolving them with a Gaussian profile to simulate the instru- 
mental broadening (FWHM=12A). A best fit using a second or- 
der polynomial gives the following result: 

logN(0 2 ) = 23.15 + 0.308 EW ma + 0.013 EW£ g70 

The zenith-corrected values of EW^io show a peak-to-peak 
range of 2.5A, with an average value of 5.46A (RMS=0.42A). 
This corresponds to log N(0 2 )=24.4 cirT 2 (2.5xl0 24 cirT 2 ), 
which is close to the value predicted by the LBLRTM model 
(24.52, or 3.3xl0 24 cm 4 ). See also the match between the 
model and real data in Fig. [9] Incidentally, as proposed by 
Stevenson (11994b . synthetic absorption spectra of this quality 
can be proficiently used to correct low resolution data for tel- 
luric features. 

To evaluate the impact of 2 variability on the extinction, 
we have run a series of LBLRTM models for different values 
of N(0 2 ) within the range deduced from the data. Although the 
variation in EW is significant (45% peak-to-peak), the effect re- 
mains confined to the two strong bands. The extinction coeffi- 
cient variation is less than 0.01 mag airmass -1 in R and /, while 
no measurable effect is seen in U, B and V passbands. The total 
contribution of 2 to the broad-band extinction is 0.01 and 0.03 
mag airmass -1 in R and /, respectively, while it is below 0.002 
in all other pass-bands^j 



10 These values have been estimated using standard passband curves, 
and assuming the spectrum of the source to be flat within the relevant 
passband. 



As far as the time evolution is concerned, the data show 
short timescale variations (cf. Sect. [63), while there is no sta- 
tistically significant signature of a long term trend over the six 
months spanned by our observations. Also, a correlation analy- 
sis between the EW of 2 and the extinction coefficient derived 
on each single spectrum in the red setting had shown that these 
are completely independent (r xy <0.2), irrespective of the wave- 
length under consideration. 

6.2.3. Water bands 

In the optical there are three bands that contribute to the extinc- 
tion, at 7200, 8200 and 9400 A, the latter being the most promi- 
nent one (Fig. [T) and affecting the z passband. The intensity 
of these absorption bands is well correlated to the Precipitable 
Water Vapor (PWV), so that a best fit using an appropriate 
radiation transfer model and suitable vertical profiles enables 
the retrieval of the PWV amount directly from observed high- 
resolution spectra (see for instance Smette, Horst & Navarrete 
12007 ). A similar method can be used to roughly estimate the 
PWV from the global EW of a given water band measured on 
low resolution spectra. This can be then translated into H 2 
column density, and finally converted into PWV. An example 
is shown in Fig. [9] which illustrates the very good match be- 
tween the data and an LBLRTM model with PWV=6mm. Using 
LBLRTM simulations run with different values of water column 
density, we have derived the following relation for the band at 
7200 A, which holds for PWV<10 mm: 

PWV ~ 0.4 EW 12 oo + 0.01 £W 7 2 200 

Before entering the EW values into this relation, one needs 
to correct to zenith the values measured at airmass X. This can 
be done through the following formula, which was derived from 
a series of LBLRTM simulations: 

EWq = EW(X) - 3.59 (X - 1) 

In principle, provided that the signal-to-noise on the contin- 
uum is sufficient (>50 per pixel), this relation can be used to 
retrieve the amount of PWV from low resolution optical spectra. 
For PWV=2mm, which is a typical value for Paranal (see be- 
low), it is £'W7 2 oo=4.1A. Unfortunately, due to the severe fring- 
ing, and the relative weakness of the feature, we could not di- 
rectly measure the amount of PWV in our spectra. Therefore, to 
estimate the variability of water features in the optical domain 
and its impact on extinction, we used the PWV data obtained on 
Paranal (Smette et al. [20071 during 2008 and 2009. The PWV 
shows a strong seasonal dependency (see Fig.[T0}, with peak val- 
ues exceeding 15mm in February-March. The distribution has a 
median value of 2.6, and in 95% of the cases is PWV<10 mm. 

Variations in the amount of PWV have a mild effect on 
the extinction coefficients, and this is limited to R and / pass- 
bands only. The synthetic broad-band coefficients deduced from 
LBLRTM simulations increase by 0.01 and 0.02 mag airmass -1 
in R and / respectively, if PWV changes from 2mm to 10mm. 

6.3. Short timescale variations 

The PARSEC data were collected with the aim of giving a sta- 
tistically robust description of Paranal extinction. Therefore, ob- 
servations have been spread as much as possible across the six 
months range covered by the project. However, on a few oc- 
casions, the same standard stars were observed at different air- 
masses during the same night, hence enabling the the calculation 
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Fig. 11. Short timescale variations observed on Dec 11, 2008 
at 5000 A. Panel A: DIMM seeing (the dots indicate the image 
quality measured on the spectra). Panel B: LOSSAM RMS fluc- 
tuations. Panel C airmass. Panel D: deviations from the average 
extinction at 5000 A. Different symbols refer to different stars: 
GD108 (circles), GD50 (triangles), and BPM16274 (squares). 
The horizontal dashed lines trace the average deviation (middle) 
and +1 sigma levels (upper and lower). 

of the specific night extinction, and the study of its variability on 
the scales of tens of minutes. The best data set was obtained 
on Dec 11, 2008, and it contains a number of repeated obser- 
vations of GD108, GD50, and BPM16274. The data, spanning 
almost the whole duration of the night, are presented in Fig. 1 1 
for ,1=5000 A. 

The extinction (on average 0.006 mag airmass -1 higher than 
the best fit value 0. 157±0.003), shows an RMS variation of 0.004 
mag airmass 



~ 1 which is comparable to the typical measurement 



error. Judging from the DIMM intensity fluctuations (Fig. 1 1 
panel B), the night was stable, with RMS flux fluctuations at 
5000 A below 1.5% on the scale of one minute. This is indeed 
reflected on the stability of the extinction, which appears to vary 
by a similar amount. We also notice that data taken at high air- 
mass do not show a systematic bias towards larger extinctions. 

7. Discussion 

The average extinction curve of Paranal, presented here for the 
first time, appears to conform to the expectations for the site. 
However, it is interesting to compare it to the data obtained for 
other two major observatories in the Atacama desert, i.e. Cerro 
Tololo (2200 m a.s.l.) and La Silla (2400 m a.s.l.), which are 
located within about 600 km from Paranal. The comparison is 
presented in Fig. 
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The data for Cerro Tololo are taken from 
the atmospheric extinction file ctioextinct.dat included in 
IRAF (Stone & Baldwin[[983j Baldwin & StoneMH- For L a 
Silla we used the table published in Schwarz & Melnick ( 1993 ), 
who reported data obtained by Tug (119771 1. These data are in- 
cluded in the atmoexan table of the MIDAS distribution. 



Fig. 12. Comparison between the extinction curves of Paranal 
(this work) La Silla (Tug [T977] ), and Cerro Tololo (Stone 
& Baldwin 119831 1. For comparison, the dashed curve traces a 
LBLRTM model computed for La Silla without aerosols. 



While Paranal and Cerro Tololo display similar behaviors 
(but see the discussion below), La Silla shows a systemati- 
cally lower extinction with the maximum deviation (~0.05 mag 
airmass -1 ) occurring at about 4000 A. Interestingly, the La Silla 
extinction curve is very well matched by an LBLRTM simula- 



tion with no aerosols (Fig. 12 dotted line). These exceptionally 
low values were noted already by Tug, who had derived a maxi- 
mum aerosol contribution of 0.01-0.02 mag airmass -1 (see also 
Minniti et al. 1989 for a comparison with two Argentinian sites). 
Since the extinction curve for La Silla was derived using data 
collected between 1974 and 1976, one possible explanation is 
that the transparency conditions have degraded in the Atacama 
desert in the last 35 years. 

Although it must be noticed that the extinction curve of La 
Silla was obtained selecting only the best nights (Tug |1977j l, the 
increased opacity is most likely due to a systematic change in the 
atmosphere above these sites. Already Stone & Baldwin (1983) 
had noticed this fact, and wrote that these observations indicate 
a non-grey increase of the extinction at Cerro Tololo, compared 
to values in use half a decade ago, averaging about 0.04 mag per 
airmass in the red and about 0.08 mag per airmass in the UV. In 
this respect it is important to note that the observations of Stone 
& Baldwin were most likely carried out before the volcanic erup- 
tion of El Chichon (March, April 1982^ which severely af- 
fected the transparency at CTIO and La Silla (see Burki et al. 
1995 for a comprehensive analysis). 

But the most intriguing aspect shown by this comparison is 
the fact that while Paranal's curve is very similar to the one of 
CTIO above 6500 A, it gradually tends to the La Silla curve in 
the blue, so that they attain very similar values below 3800 A 



11 The authors do not report the actual dates of the observations, but 
their paper was submitted in June, 1982. 
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Fig. 13. Comparison between the aerosol extinction in La Silla (1974-1976: Tug 1977'), Cerro Tololo (1980: Gutierrez-Moreno et 
al. [T982l 1982: Stone & Baldwin [T9831 1984: Gutierrez-Moreno et al. [19861 2004: Stritzinger et al. [20031 ), and Cerro Paranal (this 
work). The dashed line traces the k aer = &o/l~ a best fit law to the Cerro Tololo 1982 data (fco=0.025 mag airmass -1 , a=-1.4). 



(Fig. [12 
tion of 



The simplest explanation for this is a temporal evolu- 
aerosol mixture above the Atacama desert. 



the 

To test this hypothesis, we collected a number of published 
spectroscopic extinction curves obtained on different epochs be- 
tween 1974 and 2009. Then we computed the pure aerosol con- 
tribution as we did in Sect. [5] The result is presented in Fig. [13] 
As expected, the 1974—1976 La Silla extinction curve (Tug 
119771 ) shows an aerosol contribution which is smaller than 0.01 
mag airmass 1 . Then, the data obtained at Cerro Tololo in July- 



August 1980 (Gutierrez-Moreno et al. 119821 ) show a clear in- 
crease in the aerosol extinction, especially in the blue domain, 
where it reaches 0.05 mag airmass -1 . This increase becomes 
more marked in the next data set, obtained at Cerro Tololo in 
1982 (Stone & Baldwin 119831 ), as we said presumably before 
the El Chichon eruption. The aerosol extinction is ~0.07 mag 
airmass -1 at 5000 A, and exceeds 0.10 mag airmass -1 below 
4000 A. The wavelength dependency is well reproduced by a 
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power law with ko =0.025 mag airmass and a — -1.4 (see 
Sect. [5}, as is typical of atmospheric haze (Burki et al. |1995t . 

The next curve was obtained averaging the four data sets ob- 
tained at Cerro Tololo on January, April, May and July 1984 
(Gutierrez -Moreno et al. 1 19861 1, i.e. about 2 years after the El 
Chichon eruption. Although the data above 6000 A are a bit 
noisy, it is clear that while the extinction in the red remained 
at levels comparable to those measured in 1982, it decreased in 
the blue, making the wavelength dependency significantly shal- 
lower. This is usually interpreted as due to the evolution of the 
average particle size to larger values (see for instance Schuster 
& Parrao|200T]). 

The subsequent major volcanic event was the eruption of 
Pinatubo, which took place in June 1991. This produced very 
clear consequences on the extinction (Grothues & Gochermann 
[19921 B urki e t al. [19951 Schuster & Parrao I200T1 Parrao & 
Schuster 2003). Since then, a number of intermediate size events 
took place, none reaching the level of the Pinatubo event. The 
most relevant for our study is the explosive eruption of Chaiten 
(Southern Chile, May 2008). However, at variance with those of 
El Chichon and Pinatubo, Chaiten produced a very small amount 
of SO2, which is the main contributor to aerosol pollution. For 
this reason we tend to exclude this eruption had major effects 
on the atmosphere above the Chilean astronomical sites. In this 
respect we note that an enhancement in the SO2 content can be 
also produced by local volcanoes that are not erupting, but many 
of which have strong vapor vents. These may have variations in 
intensity of activity that are not regularly monitored or observed. 

The most recent extinction curve for the Atacama desert has 
been published by Stritzinger et al. ( 120051 1, and it was obtained 
from data taken at CTIO on seven nights in February 1999. The 
resulting aerosol term shows a rather complex behavior, with 
a rapid decrease above ~6400 A, a rather flat plateau between 
4000 and 6400 A, and a drop followed by a rapid increase be- 
low 4000 A (Fig. 13 1. Although some of the wiggles might be 



generated by the procedure used to derive the curve, the overall 
wavelength dependency is certainly different from the one ex- 
pected for a typical atmospheric haze. 

The deviation from an ordinary aerosol law is confirmed by 
our data set, which follows rather well (although with a shift 
of ~0.01 mag airmass -1 ) the points obtained at CTIO in 1984 
(Gutierrez -Moreno et al. 1986), with one remarkable difference, 
i.e. the downturn below 4000 A. A similar drop has been ob- 
served by Schuster & Parrao (2001) in their data obtained at S. 
Pedro Martir about one year after the Pinatubo eruption. Small or 
even negative power indexes have also been reported by Sterken 
& Manfroid (fT992l and Burki et al. ( fT995t as typical of vol- 
canic pollutants. We finally note that the curve by Stritzinger et 
al. (f2005 ) also shows a downturn in the blue, very similar to that 
displayed by the Paranal data. Parrao & Schuster (2003 ) have in- 
terpreted this UV drop as a possible evidence for a masking (or 
decrease) of the normal atmospheric opacity produced by the 
presence of volcanic aerosols (see also the discussion in Sterken 
& Manfroid 1 1992b . In addition, this might also be the signature 
of a more complex, possibly bi-modal, extinction law for the 
volcanic contribution (Parrao & Schuster l2003l) . 

What is surprising about these results is that such anoma- 
lous wavelength dependencies are found almost twenty years af- 
ter the Pinatubo event, when one would expect the atmospheric 
transmittance to be back to normal. But, as a matter of fact, the 
exceptional values recorded in La Silla in the 1970s were never 
observed again. The data collected in La Silla between 1974 and 
1995 clearly show that in the almost ten years that separated 



El Chichon and Pinatubo eruptions, the UBV extinction never 
reached the levels preceding the first event (see Fig. 3 in Burki 
et al. 1 1995 1 and Schwartz 2005). The extinction data presented in 
this work are on average ~0.01 mag airmass -1 lower than those 
measured at CTIO about two years after the El Chichon eruption 
(Gutierrez-Moreno et al. 1986), and ~0.03 mag airmass -1 lower 
than those published by Stritzinger et al. (2005) for the same 
observatory. This might be the signal that the atmospheric trans- 
parency is slowly approaching the high values observed prior to 
the great eruptions that took place at the end of last century. 

Further measurements in the next years will be needed to 
confirm the trend observed so far. 



8. Conclusions 

In this paper we presented the best fit extinction curve for Cerro 
Paranal, obtained combining spectroscopic data collected over 
more than 40 nights between October 2008 and March 2009. 
The main results of our analysis can be summarised as follows: 

- Above 4000A the curve is well fitted by an atmospheric 
model in which the aerosol is described by a power law 
of the form ko A a . Below this wavelength the data show 
a systematic deficit in the extinction, which exceeds ~0.03 
mag airmass -1 at 3700 A. Although this needs to be investi- 
gated with further observations, it may indicate the presence 
of pollutants of volcanic origin in the atmosphere above the 
Atacama desert. 

- During the six months covered by our observations, the ex- 
tinction distribution is characterised by a semi-interquartile 
range of 0.01 mag airmass -1 above 4000 A, with peak-to- 
peak variations up to ~0. 1 mag airmass -1 in the UV. 

- While the extinction measured at the various wavelengths 
is well correlated above 4000 A, the correlation is much 
weaker below this limit. This supports the conclusion that 
the UV portion of the aerosol contribution follows an inde- 
pendent behavior. 

- During the observing campaign, the Rayleigh scattering 
component was stable to ~0.2% (RMS). The much larger 
extinction fluctuations, seen at all wavelengths and reaching 
peak values exceeding 0.1 mag airmass -1 are attributable to 
changes in the aerosol amount and (possibly) composition. 
The median value of the aerosol contribution at 4000 A is 
0.05 mag airmass -1 . 

- Ozone was found to vary by ~5% (RMS) around an average 
value of 258 DU, in good agreement with that derived from 
the OMI data (260 DU). 

- Although the O2 A and B bands were found to vary by about 
45%, the effect on the variation of broad-band R and / ex- 
tinction coefficients is smaller than 0.01 mag airmass -1 . The 
column density of O2 derived from the data is fully consis- 
tent with the LBLRTM prediction. No seasonal effect was 
detected. 

- In the last 35 years the extinction above the Atacama desert 
has shown a marked evolution, most likely due to the two 
major eruptions of El Chichon and Pinatubo. The exception- 
ally low aerosol content measured in La Silla in 1974-76 was 
never re-established, indicating that it probably takes several 
decades before the pollutants completely fall out. 

- The usage of the IRAF extinction curve ctioextinct . dat 
leads to systematic flux overestimates of more than 4% be- 
low 4000 A. Also, it tends to over-correct by ~0.03 mag 
airmass -1 at about 6000 A, which corresponds to the ozone 
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bump. On the other hand, the usage of the atmoexan MIDAS 
table produces systematic flux underestimates that amount to 
-2% at 8000 A, and exceed 5% at 4000 A. 
- The adoption of the best fit curve presented in this paper for 
the extinction correction of spectroscopic data obtained at 
Paranal under clear-sky conditions leads to an RMS uncer- 
tainty of about 1%. 
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Fig.A.l. Slit losses for Moffat PSFs with different /3 values. 
For comparison, the Gaussian case is also plotted (dotted-dashed 
line). The upper scale indicates the FWHM seeing (in arcsec) for 
a slit width A=5 arcsec. 



Appendix A: Correction for slit losses due to 
seeing 

Because of the extended wings of the point spread function 
(PSF) typically delivered by telescopes, and depending on the 
slit width A used for the observations, a fraction of the incom- 
ing flux does not enter the spectrograph. If this loss were con- 
stant, there would be no effect on the final estimate of the extinc- 
tion coefficients. However, data are obtained under different see- 
ing conditions, and seeing tends to be larger at larger airmasses. 
Therefore, this introduces a systematic effect which leads to arti- 
ficially higher extinction estimates. Obviously, these losses can- 
not be reduced by placing the slit along the parallactic angle or 
observing through an atmospheric dispersion corrector. 

In this section we describe the method we have used to cor- 
rect the instrumental magnitudes derived from our spectra. For 
this purpose, let us introduce a coordinates reference system 
placed on the focal plane of the telescope, with its origin on 
the optical axys. Then let us indicate with P(x,y) the PSF, nor- 
malised so that: 



J PiX - 

CO xj —OO 



y) dx dy = 1 



If the slit width is A, then the slit flux losses (expressed in 
magnitudes) can be computed as: 



I P(x, y) dx dy 

00 J-A/2 



(A.l) 



for any PSF profile. For our calculations we adopted the profile 
derived by Moffat d 1969k 



P(x,y) 



1 + 
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Fig. A.2. Effect of the slit losses correction on the PARSEC ex- 
tinction at 5000 A. Upper panel: original data. Lower panel: cor- 
rected data. 



where the width parameter cr is related to the FWHM as follows: 
FWHM 



2 V2 1 // 3 - 1 



The typical values of /3, deduced from observed stellar pro- 
files, range between 2.5 and 4.0 (Saglia et al. 1993). For our 
calculations we have used (3=3, which gives a very good fit to 
the VLT FORS1 d ata in the wavelength range of interest (3300- 
8000 A). Equationl - 



A.l 



can be integrated numerically for differ- 
ent values of the FWHM, and the correction Am readily derived. 
A real example is illustrated in Fig. |A.l| For a slit width of 5 arc- 
sec (as is our case), with a seeing of 1 arcsec the losses amount 
only to ~0.008 mag, but they grow to ~0.07 mag for a seeing of 
2 arcsec. 

The effect of the correction is illustrated in Fig. |A.2| which 
presents the PARSEC extinction data at 5000 A. The uncor- 
rected data show a clear dependency of k from the seeing (the 
Pearson correlation factor is ^=0.71), while the corrected val- 
ues show no correlation (r AV =6.21). Since the seeing grows as 
A -1 / 5 (Roddier 1981), slit losses are higher in the blue than in the 
red. Therefore, not only does one expect a correlation between 
the uncorrected instrumental magnitudes and the seeing, but also 
that this relation becomes steeper in the blue. This expectation 
is confirmed by our data. An inspection of the PARSEC data at 
various wavelengths shows that the application of the correction 
removes any dependency of k from seeing and airmass. As a side 
effect, it also reduces the spread of the data points. 



Appendix B: Tabulated extinction curve for Paranal 

The merged extinction curve of Paranal is presented in 
Table B. 1 To allow the usage of the data across the whole optical 
domain, we have added four points obtained from the interpola- 
tion of an LBLRTM model. The wavelengths (8500, 8675, 8850, 
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Table B.l. Best fit extinction curve for Paranal 
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0.021 
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0.003 
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0.125 


0.003 


3375 
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0.009 
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0.003 
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0.003 


3425 
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0.007 


4725 


0.185 


0.003 
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0.002 


3475 
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0.007 
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0.182 


0.003 
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0.115 


0.002 


3525 


0.526 


0.006 
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0.176 


0.003 
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0.002 
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0.006 
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0.003 
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0.002 
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0.003 
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0.003 
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0.005 


5075 


0.153 


0.003 
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0.003 
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0.002 
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0.265 


0.004 
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(*) The last four values are interpolations from the LBLRTM model. 



and 10000 A) were selected to avoid the strong water absorption 
bands. The LBLRTM simulation was run disabling the aerosol 
calculation. Then we added the aerosol contribution using the 
best fit relation fc ael - = koA a , with fc ( )=0.013 and a=-1.38 (see 
Sect. [6X2). 



